COARSENING OF VORTEX RIPPLES IN SAND 
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The coarsening of an array of vortex ripples prepared in an unstable state is discussed within the 
framework of a simple mass transfer model first introduced by K.H. Andersen et al. [Phys. Rev. 
E 63, 066308 (2001)]. Two scenarios for the selection of the final pattern are identified. When the 
initial state is homogeneous with uniform random perturbations, a unique final state is reached which 
depends only on the shape of the interaction function /(A). A potential formulation of the dynamics 
suggests that the final wavelength is determined by a Maxwell construction applied to /(A), but 
comparison with numerical simulations shows that this yields only an upper bound. In contrast, the 
evolution from a perfectly homogeneous state with a localized perturbation proceeds through the 
propagation of wavelength doubling fronts. The front speed can be predicted by standard marginal 
stability theory. In this case the final wavelength depends on the initial wavelength in a complicated 
manner which involves multiplication by factors of 2 and rational ratios such as 4/3. 
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Uns iiberfullts. Wir ordnens. Es zerfallt. 
Wir ordnens wieder und zerfallen selbst. 
Rainer Maria Rilkej^ 



I. INTRODUCTION 
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Vortex ripples are a familiar occurrence in coastal waters, where the waves expose the sand surface at the sea 
bottom to an oscillatory flow. The name reflects the important role of the separation vortices that form on the lee 
side of the ripples in stabilizing the ripple slopes 0J^]. In laboratory experiments, one observes the formation of a 
stable periodic pattern with a ripple wavelength proportional to the amplitude a of the water motion ||[| . 

The purpose of the present contribution is to analyze a simple model that was recently introduced by Andersen 
and coworkers to describe the stability and evolution of vortex ripple patterns H . We focus here on the mathematical 
aspects of the problem, and refer the reader to the literature for further motivation and a detailed comparison to 
experiments ^ff. The model of interest is introduced in the next section. Sections 3 and 4 discuss the selection of 
the final ripple pattern for the cases of uniform and localized perturbations of an unstable initial state, while Section 
5 contains some remarks concerning the description of vortex ripples using continuum equations. 
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II. THE MASS TRANSFER MODEL 



We consider a fully developed ripple pattern of the kind obtained in a quasi one-dimensional annular geometry 
HHJ]]] (see Figure 1, which is taken from Q). Each of the N ripples is described by a single size parameter A„(t), 
which can be thought to represent its length, or the amount of sand it contains. The periodic (annular) boundary 
conditions imply that the total length L = J2n=i ^« i s conserved^. During one period of fluid motion, ripple n 
O ■ exchanges mass (or length) with its neighbors n — 1 and n+ 1. This leads to the balance equation 



^ = 2/(A„) - /(A n+1 ) - /(A„_x). (1) 

The interaction function /(A) is the central ingredient in the model. It describes the amount of sand that is transferred 
to a ripple of size A due to the vortex forming behind this ripple. 



1 From the eighth Duino elegy. 

2 Models in which the ripple lengths and masses are conserved separately have been developed in 
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FIG. 1. Experimental image of a ripple pattern obtained in an annular container under oscillatory driving. The line above 
the sand surface shows that the pattern can be fitted to an array of triangles with constant slope. The amplitude of fluid 
motion a is indicated. 

The following argument suggests that /(A) should be a nonmonotonic function with a maximum near A = a [pi. 
Small ripples create a small separation vortex which is unable to erode much of the neighboring ripples, hence /(A) 
vanishes for small A. On the other hand, even for a large ripple the size of the vortex cannot be much larger than the 
amplitude a of the fluid motion. If A 3> a, the vortex does not reach beyond the trough to the next ripple, so /(A) 
vanishes also. The mass transfer is most efficient, and / is maximal, when A ~ a. The interaction function can be 
measured in fluid dynamical simulations [|| and experiments which confirm these qualitative considerations. In 
the following we shall take /(A) to be an arbitrary, single humped function which vanishes^] at A = and at A = A max , 
and displays a maximum at A = A c . The physical meaning of A max is that new ripples are created in the troughs 
when the ripple spacing exceeds A max . Since we are concerned here with the coarsening of ripples, this process plays 
no role. 

Any homogeneous pattern with A„ = A is a stationary solution of ([!]). To investigate its stability, we impose a 
small perturbation, A„ = A + e n , and linearize ([!]) in e n . We find solutions of the form e„ ~ exp[iqn + uj(q)i], where 
the growth rate of a perturbation of wavenumber q is given by 

w(q)=2f'(\)(l-cosq). (2) 

This implies that (i) a homogeneous state is stable iff /'(A) < 0, and (ii) an unstable state decays predominantly 
through perturbations of wavenumber q = ir, in which every second ripple grows and every second ripple shrinks. The 
model (|l]) thus predicts an entire band of stable homogeneous states with wavelengths A c < A < A max , in qualitative 
accordance with experiments The existence of a multitude of linearly stable, stationary states is a property 

that the model shares with many other systems in granular physics. 

When investigating the dynamics of unstable states, (Q) has to be supplied by a rule which decides what should 
happen when the length of a shrinking ripple reaches zero. We impose simply that such a ripple is eliminated, and 
the remaining ripples are relabeled so that the earlier neighbors of the lost ripple now are next to each other. Since 
the time derivatives of the lengths of these neighboring ripples jump at the instant of disappearance of the lost ripple, 
this elimination procedure introduces a distinctly non-smooth element into the dynamics. 



III. WAVELENGTH SELECTION FROM GENERIC UNSTABLE STATES 



In this section we are concerned with the evolution out of generic unstable states. In practice this means that the 
initial state is of the form A„ = A + 8 n , where A < A c and the 5 n are random numbers smaller than A c — A. The 
evolution is then followed numerically to the point where all surviving ripples are in the stable regime, A„ > A c . 
Beyond this time no further ripples are eliminated, and hence the mean ripple length (A„) = L/N no longer changes^. 

Simulations using a wide range of interaction functions strongly indicate that the evolution selects a unique equilib- 
rium wavelength A eq which is independent of the initial conditions, and depends only on the shape of the interaction 
function. In view of the large number of linearly stable stationary states, and the deterministic character of the 
dynamics, this is a highly nontrivial property of the model, which requires explanation. While a full understanding is 
still lacking, we present here a partial solution which appears to yield at least an upper bound on A cq . 

The key observation is that (lj) can be cast into the form of an overdamped mechanical system by introducing the 
positions of the ripple troughs x n (t) as basic variables, such that A„ = x n+ \ — x n . Then the dynamics (Q) becomes 

— = (3) 

dt dx n ' 



3 Note that the dynamics (|l|) is invariant under shifts / — > /+ const. In general, A max is therefore determined by the condition 
/(A max ) = /(0). 

4 When all ripples are in the stable regime, (JILT) describes a diffusive evolution towards a completely homogeneous state. 



2 



where V is a sum of repulsive pair potentials acting between the troughs, 



V 



n+l 



E 



dA /(A). 



(4) 



This is supplemented by the elimination rule, which corresponds in the particle picture to a coalescence process^]. 

Equation (|^) suggests that the dynamics is driven towards minimizing V . It is then natural to surmise that the 
final state is determined by the minimum of V under the constraint of fixed total length L. Since the final state is 
clearly homogeneous, the quantity to be minimized is 



N I dX' f(X') 
la 



V hom (X) 

This leads to the prediction that A oq = A* , where A* is the solution of 



dX' /(A'). 



A* 



dX /(A) = A*/(A*). 



(5) 



(0) 



This is of course the analytic form of the Maxwell construction, as it would be applied to the chemical potential in a 
coexistence region. 

The prediction (j^ has several desirable properties. First, it is manifestly independent of the initial condition. 
Second, it guarantess that A* is located in the stable region, i.e. A c < A* < A max . Third, it is invariant under 
multiplication of / by an arbitrary factor^]. Nevertheless it is wrong: Comparison with simulations shows that A* > A oq 
always. Table 1 contains some typical results obtained using a family of piecewise linear interaction functions, 



/(A) 



A/A c 

(A m ax — A)/ (A max — A c ) 



A < A c 

X>X r 



In this case 



A*/A c — y A max / A c 



(7) 



(8) 



The two wavelengths A* and A eq appear to become equal, in the sense that (A* — A c )/(A cq — A c ) — > 1, when the stable 
branch of the transfer function becomes very steep, i.e. when A max — A c -C A c . On the other hand, when the stable 
branch is shallow, the prediction A oq = A* fails completely. This can be seen by considering the extreme case of an 
interaction function that remains constant for A > A c , i.e., (Q) with A max /A c = oo. Then we find numerically that the 
final wavelength is A eq ~ 1.61 A c , while clearly A* = oo. 

TABLE I. Comparison of the wavelength A* predicted by ^) with the equilibrium wavelength A cq obtained in numerical 
simulations of the model. The simulation results were averaged over 100 runs using 1000 initial ripples with lengths uniformly 
distributed in [0.5,1]. 
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An obvious interpretation of the finding A cq < A* is that the deterministic, overdamped dynamics (J3J) gets stuck in 
a metastable state before reaching the configuration of minimal "energy" V. This suggests that it should be possible 
to increase the final wavelength by either making the dynamics less damped, or by introducing noise. The first 
modification implies that the trough "particles" are supplied with a mass and a momentum variable along the lines of 
||. Preliminary simulations show that this does indeed increase the final wavelength, but not sufficiently to reach A*. 
A noisy version of the model has been described in || . Noise also increases the final wavelength, however in addition 
it introduces a new coarsening mechanism involving rare fluctuations ||, which in principle drives the wavelength 
towards A max , as long as ripple creation is not included. 



5 Related particle systems have been considered in ||. 

6 This is required because such a multiplication only affects the time scale of evolution and should not change the final state. 
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IV. FRONT PROPAGATION 



We now consider a perfectly ordered, homogeneous, unstable initial condition, A„ = A^ < A c for all n, which 
is destabilized by a local perturbation, e.g. by making a single ripple shorter or longer. Then two fronts emanate 
from the perturbed region which propagate into the unstable state and leave in their wake a stable homogeneous 
configuration at a new wavelength X^' > A c (Figure 2). The elimination of ripples occurs in the vicinity of the fronts. 
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FIG. 2. Front propagation for a parabolic interaction function /(A) = 2A — A 2 with initial wavelength A' ! ' = 1/2 and final 
wavelength A (/) = 1. The fi gure shows the ripple wavelength A n as a function of the scaled ripple number n/N (note that N 
decreases with time). 

Since the period-2 mode (q = n) is the most unstable according to (g), its growth controls the propagation of 
the fronts. Following the standard theory of front propagation into unstable states JTp| , we write the propagating 
perturbation as a traveling wave with an exponential tail, 

e„(t) = (-l) n exp[-a(n-ci)], (9) 

where c is the propagation speed and a the decay constant. Inserting this into the linearization of ([!]) we find the 
relation 

da) = J { ' (I + cosh(a)) (10) 
a 

between c and a. Localized initial conditions usually propagate at the "marginal stability" speed c* corresponding to 
the minimum of (|^|) ||lo| , and hence we expect that the front velocity is given by c* w 4.4668 /'(AW). This prediction 
is well confirmed by numerical simulations. 

We next turn to the relationship between the initial and final wavelengths. In stark contrast to the situation 
discussed in Section 3, here the final selected wavelength depends on the initial state in a surprisingly complex manner 
(Figure 3). The most prominent feature in the graph is a straight line of slope 2 which extends from A c /2 = 0.5 to 
Amax/2 = 0.75. In this regime the wavelength selection process is very simple. The growth of the period-2 mode near 
the front implies that every second ripple is eliminated, hence the wavelength doubles. When A*- 1 ' > A c /2, this is 
sufficient to bring the ripples into the stable band. We therefore conclude that 

\W = 2A« for A c /2 < A« < A max /2. (11) 
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FIG. 3. Final wavelength as a function of initial wavelength for the piecewise linear interaction function (Q) with 
Amax/Ac = 3/2. For each value of A^ a system of initially 200 ripples was simulated until all ripples reached the stable 
regime. All wavelengths are measured in units of A c 

When A*- 1 -* < A c /2, the ripples are still unstable after the doubling of the wavelength. The simplest scenario for 
the further evolution is that the new state again becomes unstable with respect to the period-2 mode, so that the 
wavelength doubles once more. Indeed a segment with slope 4 can be detected in Figure 3, which starts near A c /4. For 
smaller initial wavelengths this scenario breaks down because the accumulation of exponentially growing perturbations 
prevents the intermediate homogeneous states to become established. We have not attempted any further analysis of 
the complicated behavior seen in Figure 3 for A« < A c /4. 

A different kind of complication arises when AW > A max /2. In this case the growth of the period-2 mode terminates 
before the smaller ripples have reached zero length, because the system gets temporarily trapped in a stationary 
period-2 state with alternating ripple lengths A^ > A c and Xb < A c (Figure 4). It is possible to prove that such a 
state, which has to satisfy the constraints 

f{X A ) = f(X B ), (A a + A b )/2 = A« (12) 

always exists when A max < 2A C and A^ > A max /2. Indeed, consider the function 

F(X A ) = f(X A ) - /(A B ) = f(X A ) - /(2A» - X A ). (13) 

This function vanishes at X A — Xb — X^\ where its slope F'(X^) = 2f'(X^) > 0. Furthermore F is odd under 
reflection around A^ = A' 8 ', and F(A max ) = — /(2Aw — A max ) < 0, because 2A^ — A max > 0. It then follows by 
continuity that F has to possess two additional zeros, corresponding to a solution of (|l^) with X A > A^. It also 
follows that 

f'(X A ) + /'(A B ) < (14) 
which is the condition for stability within the space of period-2 configurations. 
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Lmax = 1.5, initial wavelength 0.9, 500 domains, t=50 
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FIG. 4. Front propagation for > A ma x/2. The figure shows a system of initial 500 ripples at time t = 50. The interaction 
function was piecewise linear with A max /A c = 3/2. Note the period-2 state appearing between the front and the homogeneous 
final state. Here A w = 0.9 A c and A (/) = 1.2 A c = 4/3A (i) . 

A stability analysis of the stationary period-2 state yields the linear growth rate 

u{q) = f(X A ) + f'(X B ) + V(/'(Aa) - /'(As)) 2 + 4/'(A A )/'(A B ) cos 2 q. (15) 

Since /'(A^)/'(A.b) < 0, the growth rate is maximal at q — tt/2, and it vanishes at q — and q = tt. We conclude 
that the stationary period-2 solution is most unstable with respect to perturbations of period 4. In effect, this implies 
that one out of four ripples is eliminated, and hence A^- 1 / A^ = 4/3. This explains the region of slope 4/3 in Figure 
3 starting around X^ « 0.8. Other rational ratios can (and do) appear in a similar manner. 



V. CONTINUUM EQUATIONS FOR VORTEX RIPPLES? 

The model ([!]) was proposed to describe the stability and evolution of fully developed ripple patterns, but it does not 
address the question of how these patterns emerge from the flat bed. In part, this reflects the fact that the separation 
vortices appear only once the pattern has reached a certain amplitude, so a different mechanism must control the 
initial instability pT| ]. On the other hand, a theoretical description that encompasses the transient evolution from 
the flat bed as well as the fully developed ripple pattern would be highly desirable, in particular for the analysis of 
two-dimensional systems Q. In this section we suggest that such a comprehensive description may be difficult to 
achieve. 

For the related problem of wind-driven (aeolian) sand ripples, a description in terms of partial differential equations 
for the (one-dimensional) continuous profile h(x,t) of the sand surface has been developed [Q. Let us collect the 
properties that such an equation should have for the case of vortex ripples under water, (i) Since the pattern does 
not depend on the thickness of the water layer, the dynamics should be invariant under constant shifts of the height, 
h — > h + const, (ii) The oscillatory driving implies symmetry under x — > —x. (iii) The slope of the ripples should 
saturate around the angle of repose, and (iv) the pattern should not be invariant under h — > —h (closer inspection 
of profiles like that in Figure 1 show that the peaks are cusp- like while the troughs are rounded, see ||). Restricting 
ourselves to terms which are polynomial in the derivatives of h, the simplest equation satisfying these requirements is 

fk = -h xx - h xxxx + (h x ) x - b(h x ) 2 xx , (16) 

where subscripts refer to partial derivatives and b is a positive constant. It is easy to see that the flat bed solution of 
( |l6| ) is unstable, with the fastest growing mode (corresponding to the initial pattern) occurring at wavelength 2ir^/2. 
The third term on the right hand side leads to a selected slope of ±1, while the last term sharpens the peaks and 
rounds off the troughs of the ripples. 

A detailed study of ( |l6| ) has been carried out by Politi Q , who shows that the wavelength of the pattern coarsens 
indefinitely as hit. This conclusion appears to apply generally to equations of the same general form with polynomial 
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terms ]12[ . Patterns which do not coarsen can be obtained only at the expense of introducing unbounded growth of 
the slope, and hence of the amplitude, of the pattern 14 1. A class of equations which contains both types of behavior 



h x 



(1 + ^)3/2 



(17) 



which arises in the context of meandering instabilities of stepped crystal surfaces |]15| . The exponent v is characteristic 
of the relaxation mechanism of the steps, the cases of immediate physical relevance corresponding to v = 1 and v = 1/2 
The analysis of this equation shows that unbounded amplitude growth occurs for — 1/2 < v < 3/2, and coarsening 
for v < -1/2. 

We therefore conjecture that local height equations generally cannot describe the emergence and evolution of 
patterns of constant wavelength and amplitude. A general proof, or the discovery of a counterexample, would be 
of considerable interest. Meanwhile, we believe that models like ([!]) can play a useful part in the analysis of such 
patterns. 
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